!*********************************************************************
! This MODULE contains a subroutine to calculate the sine and cosine
! integrals, i.e.:
!
!        x                         inf
! Si(x)= S dt sin(t)/t      Ci(x)= -S dt cos(t)/t
!        0                          x
!
! The reason it is a subroutine is that it is constructed to calculate
! Si and Ci simultaneously, as they are often both needed.
!
! Ci(-x) (x>0) does NOT return the complex part (-iPI), as this is
! seldom needed, though sici will generate a warning.
!
! Ci(0) returns a large negative number, and generates a warning.
!*********************************************************************
module mSineCosineIntegrals

        implicit none

        private
        public :: sici

contains

        subroutine sici(x, si, ci)

                use mConstants

                real(D), intent(in) :: x
                real(D), intent(out) :: si, ci
                real(D), parameter :: MAXIT = 100, & ! Max # of iterations
                    ! When x>CROSSOVER, evaluate by continued fractions method,
                    ! or when x<=CROSSOVER, use alternating series.
                    CROSSOVER = 3.0_D, &
                    ! Euler's constant, courtesy of Numerical Recipes in
                    ! Fortran 90:
                    EULER = 0.5772156649015328606065120900824024310422_D, &
!                    EPS = epsilon(1.0_D) ! Convergence criterion
                    EPS = epsilon(1.0_R4) ! Convergence criterion
!                    EPS = 1e-2_D ! Convergence criterion

                real(D) :: t, a
                complex(D) :: b, cc, bigD, f, Delta
                real(D) :: sumc, sums, sum, sign, prev, term, err
                integer :: j
                logical :: odd

                t = abs(x)
                x_is_zero: if(t == 0) then
                        si = 0.0_D
                        ci = -huge(1.0_D)

                        write(0, *) 'WARNING: Attempted to calculate Ci(0)'
                        write(0, *) '         Returning a large negative number'
                        return

                end if x_is_zero

                method: if(t > CROSSOVER) then ! Use continued fractions by
                                               ! modified Lentz's method.
                        b = cmplx(1.0_d, t, D)
                        cc = 1.0_D / tiny(1.0_D)
                        bigD = 1.0_D / b
                        f = bigD
                        j = 1
                        Delta = cmplx(0.0_D, 0.0_D, D)
                        iter_cf: do while( &
                            abs( Delta - cmplx(1.0_D, 0.0_D, D) ) > EPS)
                                maxIters_cf: if(j == MAXIT) then
                                        call warnFractionsFailed()
                                        exit
                                end if maxIters_cf
                                j = j + 1
                                a = j - 1
                                a = -a * a
                                b = b + 2.0_D
                                bigD = 1.0_D / (a * bigD + b)
                                cc = b + a / cc
                                Delta = cc * bigD
                                f = f * Delta
                        end do iter_cf
                        f = f * cmplx(cos(t), -sin(t), D)
                        ci = -real(f, D)
                        si = PI * 0.5_D + aimag(f)
                else method ! Evaluate both series simultaneously using
                            ! alternating series.
                        ! Avoid failure of convergence test due to underflow:
                        underflow: if(t < sqrt( tiny(1.0_D) )) then
                                sumc = 0.0_D
                                sums = t
                        else underflow
                                sum = 0.0_D
                                sums = 0.0_D
                                sumc = 0.0_D
                                sign = 1.0_D
                                prev = 1.0_D
                                odd = .TRUE.
                                err = 1.0_D
                                j = 1
                                iter_series: do while(err > EPS)
                                        maxIters_series: if(j == MAXIT) then
                                                call warnSeriesFailed()
                                                exit
                                        end if maxIters_series
                                        prev = prev * t / j
                                        term = prev / j
                                        sum = sum + sign * term
                                        err = term / abs(sum)
                                        series: if(odd) then
                                                sign = -sign
                                                sums = sum
                                                sum = sumc
                                        else series
                                                sumc = sum
                                                sum = sums
                                        end if series
                                        j = j + 1
                                        odd = .not. odd
                                end do iter_series
                        end if underflow
                        si = sums
                        ci = sumc + log(t) + EULER
                end if method
                x_negative: if(x < 0) then
                        si = -si

                        write(0, *) 'WARNING: Ci(x) received negative x and'
                        write(0, *) '         I''m dropping the imaginary part'

                end if x_negative
                return
        end subroutine sici

        subroutine warnFractionsFailed()
                write(0, *) 'WARNING: Sine and cosine integrals failed to'
                write(0, *) '         converge using continued fractions.'
        end subroutine warnFractionsFailed

        subroutine warnSeriesFailed()
                write(0, *) 'WARNING: Sine and cosine integrals failed to'
                write(0, *) '         converge using series expansion.'
        end subroutine warnSeriesFailed

end module mSineCosineIntegrals
